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ABSTRACT 

A simple and reliable method to study the translational relaxation of 'hot' H atoms following 
their production by chemical mechanisms is proposed. The problem is relevant to ISM, shocks, 
photospheres, atmospheric entry problems. It is shown that the thermalization of H atoms can be 
conveniently studied by a simple Monte Carlo method including the thermal distribution of background 
molecules and set the basis for further investigations. The transport cross section is determined by the 
inversion of transport data. The collisions density of H atoms in H 2 gas is calculated and discussed also 
in the context of simple theories. The application of the results to astrophysical problems is outlined 
including numerical results for the reaction H + H2O — >• H2 + OH. A simple analytical formula for 
the reaction probability during H atom thermalization is proposed. 

Subject headings: ISM: atoms — Methods: numerical — Molecular processes — Plasmas 



1. INTRODUCTION 

In many important objects atomic hydrogen H is pro- 
duced by dissociation processes from diatomic hydrogen 
H2 or other hydrogen-containing molecules like H2O, H2S 
or NH 3 . 

Exam ples are comet comas (|Id 1983T ) . Jupiter's atmo- 
sphere ()Lodders &: Feglev 201 lh and photo dissociation 
regions (PDR) (Hollcn back fc Tielens 19991 ) but the di- 
rect dissociation of H2 by electron impact is important 
also in shocks as an effect of preheating electrons, and 
in more general terms, in all non equilibrium plasmas 
produced in H2 where the electron temperature is much 
higher than the gas temperature, a situation than can 
be produced by electric fields and is easily reproduced in 
discharge plasmas. 

Atoms produced from dissociation reactions are often 
superthermal, or hot atoms, since they gain energy by 
sliding down a repulsive adiabatic potential to the disso- 
ciation limit. Hot atoms are also produced from symmet- 
ric charge exchange reactions, i.e. ^ ast + H — > Rf as t + 
H + of H atoms with fast ions of different origins. The av- 
erage energy of such atoms is a few eV. The chemistry 
of such H atoms is therefore a non thermal chemistry, 
since the usual assumption of a Maxwell-Boltzmann dis- 
tribution of the kinetic energy of H atoms cannot be held 
to relate the average thermal energy and the rate coeffi- 
cients of collisional processes. 

The production and reactions of superthermal H atoms 
is a topic of great astrophysical relevance. Many chem- 
ical reactions involving H atoms and important for in- 
terstellar chemical networks have a threshold which can 
be overcome by the availability of extra kinetic energy. 
Exampl es are the trapping of H atoms in solid hy- 
drogen (jMivazaki et al. 19831 ). the hydrogen abstraction 



from several molecules like H2O, H 2 S, hydrocarbons and 
their deuterated vers ions (|Nicholas fc Vaghjiani 1989t 
iMalcolme-Lawes 19721) . the impact ionization of I group 
metals ( Crooks et al. 19 78). It is therefore important to 
develop practical and reliable methods to characterize 
the chemical reactivity of H atoms after their production 
including the effects due to their moderation, thermaliza- 
tion and loss in impact and reactions with other species. 

This problem can be tackled by using a contin- 
uum slo wing down ap proximation similar to the Fermi 
theory (|Robson 20061 ). but a much better method is 
to apply the rigorous thermalization theory developed 
for t he case of neutron kinet ics in gaseous modera- 
tors (Weinberg fc Wiener 19581 ). 

While the resulting integral equations are very com- 
plex, a simpler approach to obtain rigorous numerical 
results is provided by the application of a Monte Carlo 
(MC) method that takes into account the thermal distri- 
bution of target particles in the collision kernel and the 
energy dependence of the cross section. 

In this paper some calculations of the translational re- 
laxation and chemical reactivity of superthermal H in H2 
are presented, with a cross section determined from the 
inversion of transport data. An accurate MC model for 
thermal particle transport developed in the past in our 
group and already validated is applied to this problem. 
The appropriate parameter set is individuated. Results 
are discussed in the light of several important astrophys- 
ical reactions. 

2. METHOD OF CALCULATION 

Calculations are performed using the MC method for 
particle transport in a thermal background described by 
Longo & Diomede (2009) and by Panarese et al. (2011). 

This method has been recently validated by com- 
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paring the calculated values of binary diffusion coef- 
ficients in different gases with calculations based on 
the Chapman-Enskog deve lopment extended to high or- 
ders ([Panarese et al. 2011h . 

Although the method is described in the above refer- 
ences, here a self-consistent short description is provided. 
The starting point is the expression of the real collision 
frequency for a H particle moving with velocity v, given 
by 



v(v) = d 3 wa (g) /(w) 



(1) 



where / is the velocity distribution function of target 
particles and a is the collision pair frequency defined as 



a g n b . 



(2) 



a is the total cross section, g = |v — w| is the relative 
speed of the collision pair and is the target particle 
density. 

The method is based on the preliminary selection of 
a maximum value for the product gcr(g) denoted by 
(gv(g))max- By replacing ga(g) with {ga(g)) max in the 
integral expression (flj this last can be rewritten into the 
form v(y) = a m ax where 



max Tib- 



(3) 



This replacement implies a potentially non-physical in- 
crease of the collision frequency, which can be compen- 
sated by using the concept of null-collision, i.e. the in- 
clusion of artificial scattering events which accounts for 
the difference a max — a but has no effect on the motion 
of H atoms. 

This solution allows an exact simple treatment of col- 
lisions in a Test Particle Monte Carlo (TPMC) model. 
This numerical method describes the motion of test par- 
ticles diluted in a bulk medium of target particles. 

In our case, the system is constituted by test particles 
of H moving in a H2 uniform bulk, in equilibrium at 
temperature T g and pressure p. 

Initially test particles are put in the origin of a three- 
dimensional space and are let to diffuse across the bulk. 
Test particles are initialized with the same energy and 
interact with bulk particles by means of binary collisions. 

For each collision, the bulk particle velocity is selected 
according to the Maxwell-Boltzmann distribution at the 
temperature T g , using a direct method of sampling. For 
this purpose, setting Vi = r sin ■& as the velocity compo- 
nent along the i-direction, a pair of values of r and t? is 
sampled from t? = 27r^i and r = (— 2kT g \n 112/ rribuik) 1 ^ 2 , 
using two random numbers 771 and 772 uniformly dis- 
tributed between and 1. Finally the value of the i- 
component of the thermal velocity in the equilibrium 
bulk is sampled as Vi — rcosT?. 

In order to remove the extra collision events used to 
equalize the collision frequency to a max , a further ran- 
dom number 773 is compared to the fraction of real colli- 
sions given by &/a max . If 773 is smaller than this quantity, 
the collision is effective. 

After an effective collision, the relative velocity vector 
must be rotated according to two polar angles, namely 

the scattering angle, and <p, the azimuthal angle. This 
last is uniformly sampled in the interval [0,27r], while the 
selection of $ depends on the interaction model. 



Once the scattering angle is known, the scattering is 
treated taking into account the correlation with the old 
particle velocity using Euler angles: the relative velocity 
vector after the collision, g* , is calculated as 

g* = g x cos d + B sin tp sin t?, 
g* y = g y cos -d ~ B~ x sin -d (gg z cos tp + g x g y sin <p) , (4) 
g* = g z cos i9 + B~ x sin d (gg y cos ip - g x g z sin tp) , 



where B = {g 2 ^ 



■i V < 9z) and g = (9x,g y ,gz) as above is 
the relative velocity before the collision. 

t9 is determined from a quadrature of the interaction 
potential (j)(r) based on the known value of the impact 
parameter b 



tf(b,E) = tt-26 



1-1/2 



E 



dr 



where r m is the distance of closest approach. The value 
of b is obtained from b — b maxy fffy. <j> is simply given by 
4> = 2^774. In case of isotropic elastic scattering, cos (7?) = 
1-2775. 

The motion of the colliding particle of mass m c relative 
to the bulk target particle of mass m* is equivalent to the 
motion of a particle of mass /1 = m c m t / (m c + m t ) rela- 
tive to a centre of force. The collision energy is calculated 
as a function of the relative speed g of the interacting pair 
by E C oii = /A9 2 / 2 - 

As in a binary interaction the centre of mass velocity 
is a constant, the velocity of the colliding particle after 
the collision is given by 



77l c V + m t w 

777 c + 777 t 



rrii 



m t 



(5) 



The time difference between one collision (including 
null collisions) and the next one is given by the formula 



At 



Inrje/ar, 



(G) 



where 776 is again a random number from a uniform dis- 
tribution, < 77 6 < 1. This time is inversely propor- 
tional to the gas density. An appropriate parameter to 
measure the degree of thermalization of a H atom in the 
bulk medium is therefore provided by the product nr, 
where r is the average lifetime of a H atom before its 
chemical or diffusion loss. Calculations performed using 
different values of n and r, but leaving the value of the 
product nr unchanged, produce the same results. 

The finite lifetime r is obtained in the simulation by 
removing the injected particles with a probability 1 — 
exp(—At/r) before each collision. The removed particle 
is labelled in such a way that the computer simulation is 
not accounting for it anymore. The simulation proceeds 
until all particles are removed. 

This procedure is exact when r is not dependent on 
the atom speed and can be used when this approxima- 
tion is considered feasible. In specific cases r will be 
determined by the chemical network assumed. In molec- 
ular clouds, collisional loss reactions with and II ~ 
and other species can be important, while in the labora- 
tory H atoms are often lost in reactions with purposely 
added scavengers uncommon in space, e.g. iodine atoms. 



MC translational relaxation of H atoms 
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For these collisional losses r = (^l^l)" 1 where is 
the appropriate rate coefficient while is the number 
density of the corresponding reaction partner. The pho- 
toionization or the radiation ionization of H atoms are 
well described by a constant r as well. Heterogeneous 
recombination on reactor walls in the laboratory and on 
dust grains in space can be treated similarly, but the for- 
mula must account for diffusion times. In case r has a 
strong dependence on the H atom speed or if a loss cross 
section ctl is used instead of a loss rate coefficient, the 
removal process can be improved. For any collision a 
candidate collision partner for any loss channel is sam- 
pled from a Maxwell distribution. The probability that 
a collision results in H atom loss is pl = ogn^ /a max and 
the loss channel is selected if 777 < Pl- 

Collision events are sampled on a uniform grid on an 
axis representing the variable u = Ln(e/eo). In this way 
the quantity sampled is the collision density P(u)du. 
This quantity is normalized to one atom produced per 
unit time. 

3. RESULTS AND APPLICATIONS 

The role of hot atoms in chemical reactions is well 
known as mentioned in the introduction. With the extra 
energy provided, these atoms can increase considerably, 
sometimes by orders of magnitudes, the rate of chemi- 
cal reactions in which they are involved. There are two 
ways in which the non equilibrium translational kinetics 
of H atoms can be described, as a result of calculations 
by the method described in the previous section. One 
possibility is to calculate the so-called collision density, 
or equivalently the energy distribution. The second one 
is to calculate the rate coefficients of relevant reactions 
with trace species present in the bulk gas (H2). 

The collision density P(e) is defined in such a way 
that the number of collisions of each atom in the en- 
ergy range (e, e + de) is given by P(e)de. The knowledge 
of this function allows to calculate any collisional rate 
once the probability lUj(e) of accessing the i-th channel 
is known. Alternatively, the usual kinetic distribution 
f(v) can be calculated, if relevant, by the expression 
f(v) v(v) = P(mv 2 /2) mv where v(v) is the collision 
frequency for test particles of speed v reported at the 
beginning of the previous section. The collision density 
has been the subject of classical theoretical treatments 
(see also below) but in the case of slightly superthermal 
particles decelerated until thermal energy (thermaliza- 
tion), which is relevant here, a full numerical approach 
is more appropriate. The collision density has been the 
subject of stochastic calculations, beginning from Rebick 
& Dubrin (1970). In their work the authors provide the 
first results for H atoms in Xe, a diluent selected in view 
of the high molar mass and relative chemical inertness. 
They neglected the effect of the target speed in some 
steps of the calculations, which is a quite acceptable ap- 
proximation for their test case where the mass ratio is 
very high, but not appropriate for H in H2 . The collision 
frequency and dynamics is exactly accounted using the 
method of this paper. 

In order to calculate the frequency of H-H2 collisions 
a reliable momentum transfer cross section is necessary. 
Such cross section can be evaluated by fitting the binary 
diffusion coefficient D of H atoms in H2 gas at different 
temperatures. This last is calculated using the collision 
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Fig. 1. — Fit of the binary diffusion coefficient of H in H2 cal- 
culated using data by Stallcop et al. (1996), and using the model 
cross section in the text. 

integrals reported by Stallcop et al. (1996) and based 
on an accurate interaction potential. The MC method 
described in the work by Panarese et al. (2011) is used 
to obtain D as a function of T g from a guessed cross 
section of the form <jq(vo/v) p with vo = 3.75 km/ s. The 
best fit (fig. [J) is obtained for 00 = 17.5 A 2 and p = 0.9. 
This values are used in the subsequent calculations of H 
thermalization. 

Fig.[T]is cut at T g = 2x 10 3 K, since above this temper- 
ature, the equilibrium dissociation of H2 at latm cannot 
be neglected. However, the comparison of the two values 
of the formal quantity D(T S ) at higher T g is still mean- 
ingful as a check of our simplified cross section at high 
energy. These two values are actually very close also at 
T g ~ 10 4 K. 

The numerical parameters revelant to specify the ther- 
malization regime in our problems are: the initial H en- 
ergy eo, the gas temperature T g , the gas density n and 
the average lifetime of H atoms in the gas, r. This last 
parameter describes the effect of processes which remove 
H atoms from the gas. This process can be e.g. the pho- 
toionization of H atoms. In this and similar cases r is 
defined as usual through the photochemical expression 
t = riH /Rate where Rate is the chemical loss rate of H 
atoms. 

In order to reduce the number of parameters, in this 
study eo is fixed to 2.4 eV. This value is compatible 
with direct dissociation from the 1 T, g ground state to the 
triplet state correlated to H(ls) + H(ls) and provides an 
important example. This choice of a delta function for 
the energy distribution of the source can be extended to 
more complex cases by simple statistical sampling of the 
more involved source function. 

The collision density in case of H atoms originated with 
a single energy is characterized by a singular contribu- 
tion corresponding to the first collision of H with H2 with 
a H energy by necessity equal to eo- A finite amplitude 
for P(u) could be obtained by selecting the initial energy 
from a Gaussian like in the work by Prisant et al. (1978). 
Follo wing the custom for P(u) p lots in nuclear applica- 
tions (|Weinberg fc Wiener 19581 ) this contribution have 
not represented. 

Results for P(u) are better discussed in the light of 
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Fig. 2.— Solution in the cases of T S =1K, T 9 =20K and 
T 9 =100K. The Placzek function, characterized by a discontinu- 
ity for e=eo/9, is evident. The singular contribution at e = eo (see 
text) is not represented. 

the semi-analytical theory of Placzek (1946), based on a 
iterative solution of a simplified integral equation, which 
is obtained for P(u) in the limit of zero gas temperature 
and for rigid sphere elastic scattering. It is important 
to note that stochastic approaches allow to remove both 
limitations of Placzek equation. Nevertheless, the solu- 
tion of this equation allows to establish three important 
general features which are approximately valid also for 
T g > 0: 

1. In the so called asymptotic region far below eo but 
close to the bulk thermal energy, P(u) = l/£ where 
£ is the average logarithmic energy loss for single 
collision. For H in H2, this means P(u) = 1.4. 

2. Very close to the source energy eoj P( u ) is given by 
a simple expression as a function of the mass ratio 
A = ni H2 /m H = 2 here, i.e. P(u) = (A + l) 2 /4yl = 
1.125. 

3. A discontinuity is present in the function at e = 
e ((A - I) /{A + l)) 2 = 0.11e associated to the 
lowest energy value of the initial redistribution by 
the first collision. In fact, the maximum loss for the 
H energy corresponds to a head-on collision with an 
H2 at rest: in this case the final H speed is equal 
to (A — 1)/(A + 1) = 1/3 the initial H speed. 

Such features are confirmed by our calculations for dif- 
ferent values of T g , as shown in fig. [5] In this case the 
particles were removed when moderated below 0.024 eV 
in order to avoid any interference in the moderation ki- 
netics. This choice produces a non-exponential loss of 
the particles during moderation. 

In fig. [3] the collision density corresponding to T g = 
100A" and different r is reported. The equilibrium result 
is also reported for reference. The plots illustrate the role 
and effects of the two fundamental parameters t/tq and 
T g on the calculated collision density. r = (nero w o) _1 is 
the characteristic collision time. An estimate of To for a 
typical ISM density of 10 3 cm~ 3 leads to r ~ 1.5 x 10 6 s. 

A superthermal tail is associated to the slowing down 
of freshly produced H atoms, whose distribution shape 



is affected by T g and the shape of a g but essentially cor- 
responds to the Placzek solution for A — 2. The low 
energy component of P(u) is essentially determined by 
the equilibrium contribution, the relative importance of 
the two being controlled by the parameter t/tq. Even 
at energies much higher than those corresponding to the 
Maxwellian bulk, in the tail region, the thermal distri- 
bution of molecules cannot be neglected. This has the 
effect of smearing the Placzek peak of fig. [2] which does 
not appear anymore as a salient feature at temperatures 
higher than a few 10's K even for rigid sphere scattering. 
Only detailed calculations can then establish accurately 
the tail shape. 

The presence of the superthermal tail is expected to 
limit the validity of the traditional chemical kinetics of 
astrophysical objects and phenomena in several ways. 
For example, the rate coefficient of a reactive process 
involving H atoms will not be given by the equilib- 
rium formula based on the Maxwell-Boltzmann distri- 
bution, but a superthermal contribution to the rate co- 
efficients will ne ed to b e accounted for, i.e . K to t = 
K eq (T g ) + K neq (jlp 1983t [Crooks et al. 1971 ). This last 
contribution K neq can be determined in specific cases by 
running a MC simulation based on the prescriptions pro- 
vided in this paper. 

Details differ depending on the gas temperature and 
the threshold of the reaction. For reactions with no 
threshold (where effects are still possible due to cross 
section shape or resonances) or a threshold comparable 
to the thermal energy, the rate coefficient can be cor- 
rected by the modification of the translational distribu- 
tion which is not described by a Maxwell-Boltzmann law. 
The expression of the rate coefficient K p for a process p, 
whose related cross section is cr p , is a functional of the H 
translational distribution 




^(e)e 1/2 / cm (e)^, (7) 



where f cm (e) is the distribution of impact energy in 
the cm. frame. In the equilibrium case f cm (e) is 
a Maxwellian distribution with m replaced by \x and 
K p is a function of T g . When a superthermal com- 
ponent is present, this equation is still valid but K p 
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Fig. 3. — Collision density vs particle energy in the case T S =100K 
for various values of t/tq . The singular contribution at e = eo (see 
text) is not represented. 
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cannot be written as a function of the gas tempera- 
ture. Examples of reactions which may be reconsid- 
ered in this light a re the radiative association reac- 
tion s of H with C + flBarinovs Sz van Hem crt 200(f) and 
H + (jStancil et al. 19 93) which are of importance for the 
interstellar medium and for the early Universe chemistry 
respectively. 

Superthermal atoms can also affect the cooling func- 
tion of the gas, due to the strong energy depen- 
dence of the collisional deexcitation coefficients, i.e. 
the rate con stants of the react ions B+B 2 (u,j) — > 
H+HaCz/,/) (jCapitelli et al. 2006( ). The complex ther- 
malization of H atoms will also affect preheating in strong 
shocks in hydrogen. The atomic component is more mo- 
bile than the bulk gas because of the much higher dif- 
fusion coefficient, and can contribute to affect the shock 
profile. A similar effect has been studied in the past by 
Bruno et al. for nitrogen sho cks by a Direct Sim ulation 
Monte Carlo (DSMC) model (jBruno et al. 20021 ) . 

A much more important and complete nonthermal ef- 
fect is expected in cases where the reaction threshold is 
much higher than the average thermal energy: the ther- 
mal rate of such processes is essentially zero, and calcu- 
lations based on the method described in this paper can 
be used to define a specific rate Q which enters in an 
expression like 



Reaction rate = Qxsdnu / dt 



(8) 



where Xs is the mole fraction of the trace partner s and 
dun I dt is the production rate of high energy atoms per 
unit volume and time. 

Q is therefore the number of reactive collisions of H 
atoms with a reactive trace species during their moder- 
ation, normalized to unit fraction of the trace (but in 
practice the fraction is << 1 being s a trace). 

Examples of such reactions are the reaction of hydro- 
gen with B~ leading to detachment in the so-called non 
associative channel: 

BT+H — >2B + e. (9) 

This process has a threshold of about 0.7 eV and its 
cross section increases by an order of ma gnitude in the 
energy range 1 — 10 eV (jJanev et al. 20 03?). This reaction 
can be of relevance for chemical networks connected with 
radiation transport in astrophysical systems in view of 
the great importance of B~ ions as radiation absorbers. 

Another example is the abstraction of B atoms from 
water, e.g. the reaction 



E + B 2 



B 2 + OB. 



(10) 



This is potentially an important reaction in astrophys- 
ical systems, since it represents a source of hydroxyl rad- 
icals in rotovibrationally excited states. In spite of this, 
the above reaction is usually not included in chemical 
networks in view of its relative high threshold, about 
1 eV. This situation, however, can change in presence of 
superthermal atoms, and for future reference a calcula- 
tion of its rate coefficient is provided as a function of the 
initial energy of B atoms. 

The rate constant Q(eo, T g ) can be evaluated by calcu- 
lating, in a MC simulation, the sum 



Q(eo,T g ) 



jy-1 9H/H 2 Q Vabs(gH/H 2 o) ^ 
^ 9H/H 2 Vel{9H/H 2 ) 



where a a b s is the abstraction cross section and a e i is the 
cross section of the elastic process, N is the number of B 
atoms in the simulation and the sum includes all E/E 2 
elastic collisions. 

The cross section cr a b s (e) is that reported as a full line 
in fig. 4 in the paper by Brouard et al. (2004). 

The results in fig. S] show that hot B atoms can be 
an effective source of OB radicals in presence of traces 
of water diluted in B 2 . The slight increase observed for 
the highest value of T g (2 x 10 3 K) is due to the effect 
of the relative speed in B/B 2 collisions on the energy 
distribution of B atoms and confirms the necessity of 
calculating accurately the effects of the relative speed 
under some circumstances. 

The use of the MC method described here produces 
an exact evaluation of the quantity Q. An approximate 
analytical estimate can be provided under the hypotheses 
that: 

1. the Placzek equation applies; 

2. the reactive cross section in the energy range (0, eo) 
is larger for e close to eo; 

3. the mass of the B atom is much smaller than that 
of the reaction partner. 

The last hypothesis allows to neglect the difference be- 
tween the energy of the atom and the energy avail- 
able for the collision. An analytical solution of the 
Placzek equation in the energy range (aeo,eo), where 
a = ((1 — A) /(I + A)) 2 , is given by the expres- 
sion (jGlasstone et Edlung 1951 : 



o/(l-a) , 
P(e)=5(e-e ) + ^- 



(12) 



I- a eV(i-a)" 

Baving defined the function 

q(e) = a(e)/<r el (e), 

Q can be written as the result of the integration in the 
range (0, eo) of the product qP, i.e. 



(13) 



a/(l-a) r e 



Q = ?(eo) 



L 



1 - a 



dc 



717(1^)' ( 14 ) 



T 




equation (1 4) 

Tg=50K 
Tg = 200K 
Tg=500 K 
Tg = 2000K 



"T" 



"T" 



"T" 



T 



1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 
Fig. 4. — Q vs particle energy in the case T g = 50K, T g = 200 K, 
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where the lower integration limit is a reasoble cutoff (sec 
n g EJ), since at such low energies the collision density is 
thermal. 

In fig. E] results of eq. (fl4| for several values of eo are 
reported for the reaction of H abstraction from water. In 
this case the agreement is very good, since this reaction 
meets well the requirements 2 & 3 above. Of course the 
results at the highest value of T s are not matched since 
at high T g the Placzek equation is not exact anymore. 
Eq. (fT4]l can be used as an alternative to MC calculations 
for fast, approximate calculations in other cases. 

4. CONCLUSIONS 

In this paper it has shown that the problem of the 
thermalization of hot H atoms produced by photochem- 
ical, electron or ion impact processes in hydrogen based 
plasmas of astrophysical relevance can be conveniently 
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addressed by a simple MC procedure recently developed, 
which, while being very simple, allows for a rigorous 
treatment of the thermal distributions of background 
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